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ABSTRACT 


A computer simulation model is developed treating wave 
propagation in a random, inhomogeneous ocean as transmission 
through a linear time-invariant, space-variant random communi- 
cation channel. The ocean volume is modelled by an index of 
refraction which is decomposed into a depth-dependent deter- 
ministic part and a depth-independent Gaussian zero-mean 
random part. Computer simulated output electrical signals 
were generated that depend on the complex frequency spectrum 
of the transmitted electrical signal, the far-field beam 
pattern of the transmit array and the random transfer function 
of the ocean medium. Output was generated for different test 
cases. In all cases the transmit electrical Signal was repre- 
sented by a finite Fourier series and random cases were 
modelled by a random number generator. The computer simulated 
Output electrical signals were then processed by a 3-D DFT 
beamformer and the results for the deterministic inhomogeneous 
and random inhomogeneous cases were compared to the homogeneous 
non-random case in order to study the effects of the medium 


On signal distortion and source localization. 


CABLES CONTENTS 


i INTRODUCTION ------------------------------------ 
aol. Peony SON Siar ece PULLERS SEMULATTION MODEL =—=-------- 
A. MODEL DESCRIPTION --------------------------- 
1. Geometry of the Model ------------------- 
2. Transfer Function ----------------------- 
3. Output Electrical Signal ---------------- 
B. CONSTRAINTS ON THE MODEL -------------------- 
Jue Geometrical Constraints ----------------- 
2. Sound Speed Profile -----~---------------- 
3. Frequency Spectrum of Input Electrical 
Signal ---------------------------------- 
TII. IMPLEMENTATION OF THE MODEL --------------------- 
ae NUMERICAL METHODS --------------------------- 
ie, Piast loon. gles. ——<—<—<—_$——=——$—— — ——— ——— 
Ze Pee eartOn OL simceqrdcion Methods —=—-—-—- 
B. PROGRAM OUTLINE ----------------------------- 
JES SIMULATED DATA ---------------------------------- 
A. PRESENTATION OF COMPUTER SIMULATED DATA ----- 
B. TEST CASES ---------------------------------- 
C. RESULTS ------------------------------------- 
V. CONCLUSIONS AND RECOMMENDATIONS ----------------- 
LIST OF REFERENCES ------- ------------------------ 


Pia TEAL DISTRIBUTION LIST 


12 


2 


eZ 


16 


18 


Zo 


20 


Za 


22 


24 


24 


27 


Sk 


45 


aS 


a 


50 


62 


88 


oo 


2 


eee 


DV 


TVe3 


Iv.4 


LV <2 


6 


De 


IV a8 


ivi? 


EO 


EV je dek 


LISt oe ab oe. 
Magnitude of Frequency Response Comparing 
Cases Which Differ Only in Range ----=------<----— 69 


Magnitude of Frequency Response Comparing 
Cases Which Differ Only in Cross Range --------- 69 


Magnitude of Frequency Response Comparing 
Cases Which Differ Only in Frequency ----------- 70 


Magnitude of Frequency Response Comparing 
Homogeneous and Inhe@m@qenCoUcm Cac Sete ee 70 


Magnitude Of Yonld, vans.) ence mile 


Cosine u for HMG1l, INHMG1 ---------------------- 7 

Magnitude of Youlq,t yn) vse ite elem 

Cosine u for HMG2, INHMG2 ---------------------- Ta 

Magnitude of Yeu(q,m)S) vous Direceren 

Cosine v £0r° HMG MG] a iZ 

Magnitude of Yonlq,m,S)}) vs, sDeneerton 

Cosine v for INHMG]1, INHMG2 -------------------- 72 

Estimated Values of u_,v_,8_,W_ ---------------- 73 
0 oF 4O © 

Estimation of Fourier Series Coefficients 

You (G70,0) for HMG8, INHMG8 and RINHMG8 -------- 74 

Estimated Values of u_,v_,8_ and w_ ------------ 74 
oO 1 © O 


hoe 


‘i 


He. 


nS 4 


nA 


Ss 


6. 


1 ae 


nS 


ro 


Lister (GURES 


Model Geometry ----------------------------------- 


Oscillatory Integrand Showing Decreasing 
Contributions to the Integral -------------------- 


iyeteal Behavior Of Integrand in Eq. (3.25) 
Around Stationary Point -------------------------- 


Flowchart of Main Program Module ----------------- 
Flowchart of Module CALCHF ----------------------- 
Calculation of Subintervals by CALCHF ------------ 


Estimation of Direction of Propagation 
for Case HMG]l ------------------------------------ 


Estimation of Direction of Propagation 


for Case HMG2 ------------------------------------+ 


Estimation of Direction of Propagation 
for Case INHMG] ---------------------------------- 


Estimation of Direction of Propagation 
for Case INHMG2 ---------------------------------- 


EoeimMatLon of Direction of Propagation 
for Case RINHMG] --------------------------------- 


BStimation of va for Cases RINHMG2 and 
RINHMG1 (Different Seeds) ------------------------ 


Angular Spectrum for Cases INHMG2 and RINHMG2 ---- 


Angular Spectrum for Case INHMG8 at f = 1000 
Te Gee ee ee ee ee Se 


Pagular Spectrum for Case INHMGS at £ = 3000 Hz -- 


Angular Spectrum for Case INHMG8 at f = 4000 
and f = 5000 ee ee ee ee a ee ee a ee 


Angular Spectrum for Case RINHMG8 at f = 1000 
aneceet =| 2000 Hz8—=]{—={—-—{—{—<—{——— —- —-_—- - 2 


Angular Spectrum for Case RINHMG8 at f = 3000 
iG el Oleh = oe ee Ke 


Angular Spectrum for Case RINHMG8 at f = 5000 Hz - 


80 


eal 


82 


83 


84 


85 


86 


7 


ACKNOWLEDGMENT 


The author expresses his sincere thanks to Professor L.J. 
Ziomek. His generous aid and encouragement have been of great 


support in completing this thesis. 


ieee eR Omue tT TON 


Based on the linearized acoustic wave equation for small 
amplitude acoustic signals, sound propagation through the ocean 
medium can be modelled as a signal passing through a linear 
time-variant, space-variant random filter. The term "space- 
variant" implies that the sound speed profile is a function 
of position. The space variant property results in scatter 
Semangular spread due to refraction. If the filter is 
"Space-variant," then an isospeed medium is implied. As a 
result, there will be no refraction, and hence, no scatter or 
angular spread since the sound rays will be travelling in 
straight lines. The term "time-variant" implies motion of, 
for example, discrete point scatterers, the ocean surface, 
and the transmit and receive arrays (apertures). The time- 
variant property results in both Doppler spread and spread 
in time delay. 

By using a linear system's theory approach, different 
propagation paths, such as direct paths, bottom reflected 
paths and surface reflected paths, can be modelled as the 
outputs from linear filters. Furthermore, different trans- 
mit Signals and transmit and receive far-field directivity 
functions can easily be coupled to various transfer functions 
of the random, inhomogeneous ocean medium in order to study 
their effects on signal distortion and source localization 


using various space-time signal processing algorithms. 


Work based on this approach has been done by Laval [Refs. 
5,6], Laval and Labasque [Ref. 7], Middleton [Refs. 8,9] 
and Ziomek [Ref. 1]. Middleton [Refs. 8,9] described propaga- 
tion phenomena in a random inhomogeneous ocean with the use 
of space-time operators but did not concern himself directly 
with the derivation of random, time-variant, space-variant 
transfer functions. 

Ziomek [Ref. 2] derived a time-invariant space-variant 
random transfer function of the ocean volume based on the WKB 
approximation. Ziomek [Ref. 3] also derived an equation for 
the random, output electrical signal at each element in a 
receive planar array of complex weighted point sources in 
terms of the frequency spectrum of the transmitted electrical 
Signal, the transmit and receive arrays, and the random ocean 
medium transfer function. 

The purpose of this thesis is to implement on a computer 
the equation for the output electrical signal as derived by 
Ziomek [Ref. 3]. This mathematical computer model simulates 
the effects of wave propagation through a random inhomogeneous 
ocean. The ocean volume is modelled by the transfer function 
derived in Ziomek [Ref. 2] and incorporates an index of 
refraction which is a function of depth. The index of refrac- 
tion is decomposed into a depth-dependent deterministic com- 
ponent and a depth-independent zero mean Gaussian random 
COMPONeENe. 

Section II is devoted to a theoretical description of the 


computer model. The notation and geometry of the computer 


EG 


model are explained and the fundamental equations on which 
the computer model is based, are stated. Also, some additional 
constraints are made to allow for simpler implementation. 

Section III describes the numerical techniques used and 
the resulting mathematical procedures which are implemented 
in the program. Furthermore, a brief outline of the computer 
program iS given. 

In Section IV a number of test cases are defined which 
were used to test the computer model. Computer simulated 
output electrical signals were generated for test cases 
incorporating deterministic homogeneous, deterministic 
inhomogeneous and random inhomogeneous medium models. The 
output signals were processed by a 3-D DFT beamformer and 
the results for the deterministic homogeneous case were used 
as the baseline case. Both deterministic and random 
inhomogeneous cases were compared to this baseline case in 
order to study the effects of the medium on signal distortion 


ema source localization. 


Jel 


Il. THEORY ON THE COMPUTER SMU TlOnmiemrne 


A. MODEL DESCR Eire) 

The mathematical computer model calculates the output 
electrical signal at each element of a receive planar array 
of point sources in terms of the complex frequency spectrum 
of the transmitted electrical signal, the transmit and receive 
arrays, and the random ocean medium transfer function. The 
theory for the model is based on time-variant, space-variant 
random filters as discussed in Ziomek [Ref. 1]. JZiomek [Ref. 
2] derived a time-invariant space-variant transfer function 
for a random ocean volume where the index of refraction, which 
is composed of a deterministic and a random component, is a 
function of depth. This transfer function is based on the 
WKB approximation, which is a solution to the wave equation 
when the index of refraction is a function of depth. Further- 
more, Ziomek [Ref. 3] derived an equation for the random 
output electrical signals appearing at each element in a 
receive planar array of complex weighted point sources in 
terms of the frequency spectrum of the transmitted electrical 
Signal, the transmit and receive arrays, and the random ocean 
transfer function. The computer model combines and implements 
the equations derived by Ziomek [Refs. 2,3] for the transfer 
function and the output electrical signals. 

iy BGeeirere 7 lommene my ocer 

Figure 1 shows the geometry of the communication 


channel as simulated by the computer model. The model 
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describes the wave propagation between a transmit planar array 
and a receive planar array whose centers are located at 

(Xo VY G12) and Vee creat respectively. The x, y and z axes 
describe cross range, depth and range, respectively. The 
transmit array, which is a parallel to the XY plane, consists 
of MxN complex weighted point sources, where M and N are odd. 
The sources are equally spaced in the X- and Y-directions 

with spacings of dq. and a meters, respectively. Note that 

in general dq. “ - The complex weights are assumed to be 
separable and given by 


Cae We a, exp {63 (2.1) 


for the weighting in the X-direction and by 


is b_exptjo} 2) 


for the weighting in the Y-direction. The combined weight 
factor for each transmit element is eee where m and n are 
the indices describing the transmit element at position 

oS ao The complex weights are used for ampli- 
tude shading and beamsteering of the transmit beam pattern. 
The receive planar array, which is also parallel to the XY 
plane, consists of M' xN' elements, where M' and N' are odd. 
The sources are equally spaced in the X- and Y-directions with 
Spacings of qd and a meters, respectively, where in general 
dy # ve The position of a receiver element with indices 


(1,0) Is Given by (rn tidr ry tad) 2) - 


14 


The speed of sound in the medium is assumed to be 
only a function of depth and is given by c(y). The reference 
speed of sound Co is the speed of sound at the center of the 
meansmit array, 1.ée., ce = cly.)- The index of refraction 


n(y) is a function of depth and is composed of a determinis- 


tic component and a random component: 
nly) = ¢ /ely) = npyly) + n, Cy) (23) 


where ny is the deterministic component and Np is the random 
component of the index of refraction. The initial propagation 
vector will be denoted by Kae ie = 2mf/c. is the magnitude 

of the initial propagation vector at depth y = Ves Cise np at 
bm@e center of the transmit array). The components of the 


initial propagation vector along the x, y and Zz axes are 


Seseribed in terms of the direction cosines Uos Ms and Wo? 


k eral (24) 
X © oO 

k OR 7 5) 
y @ © 

k en ae (25.6.) 
Zz Oo Oo 

and 

we = 7 - yt (2) 

O Oo O 


Note that the magnitude of the initial propagation vector at 


the transmit elements which are located at a depth BS 4) MBE 


iS 


from the center of the transmit array is, in general, not 


egual to k, and is given by Ky? 


KA = io ase C2 c)) 


One can also use spatial frequencies (in cycles/meter) to 


denote directions of propagation: 


- = US Ai 18 = WRT rE = w/a (24938 


2. dranstes Weunetien 
The medium transfer function Hy 1s given by a functicn ye: 
frequency f, spatial frequency fy and receiver depth y for 


a given source depth V5 [REG wo Zhe 


Hy ltt zy) = Ay(f,f iylexpij lOyp (t,t iy) 


+ Sup (Ef iy) d3 (2.13 


where Ay is the amplitude term given by 


a 22 -l1/2 
Ay = ie) = ves) (221i) 
and oup? Our are the deterministic and random phase terms, 
respectively: 
= (-k*/4n£ [n2 ljd 1) 
‘mp = (Ko/4tEy) J Nyc) -1lidz C2. 
O 


Ts 


ee eaic (ze ) 


Mrpeeetransrer function is valid in the absence of turning 


points and when the following inequality is obeyed: 


In*(y) - lj/ve < 1 (2.14) 


In the case of a speed of sound profile with a constant 
Gradient g, the deterministic component of the index of 


refraction becomes: 
Nyly) = en) (crea y=y)) (2205) 


Substitution into equation (2.12) results in 


ee 
ee Sa Ss = = 
8D = oe ope bea oy ne (2.16) 
As far as the random phase term is concerned we will 
proceed by assuming that the random component of the index of 


refraction is not a function of depth and can be represented 


by a Gaussian, zero mean random variable with variance ae 
np (Z) ee oi C2217) 


where Np is nea) and Nye which is the normalized random 


component of the index of refraction, is N(0,1l). When Np 


foenOtea £UmMeEiOon Of depth, equation (2.13) becomes: 


ey 


w~ 


5 y 
Cn =) ae f Heal (2258); 


O 
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Substituting the equation of n. for a COnStant gradation eens 


D 
equation (2.18)9 teadstee- 
we So 
oe = 7 =I ne Int{n, (y) } (2 13 


For a homogeneous medium, equations (2.16) and (2.19) become: 


oD = 0 (2... 2038 


w 


je eee! a 
CA 7 ne ly -y,) (2.21) 


Equations (2.11) through (2.13) for the anolaeucera 
phase terms of the transfer function show that the méedizunma 
the effect of angle modulating the transmitted sound field, 
also referred to as Scateenind. 

oo) CUE OU ERE Hecumicc ime oma 
The output electrical signal at a receive location 


PSEC cS nil gt) 1S given by [Ref. 3]: 


00 1. a 
yltx Hidtjytadyj2.) = f xe) ff f 
—co -l a 
q 
(NeyZ 
ceny 79 Sr So a EP) DTV fg JOC, (IU Va) : 
n=- (N- 
(M-1) /2 ow 
) Cc exp{—jk ux, Idv du exp{j2nftid; ujtv, <1 (2.22) 


m=-(M-1)/2 ™ 


is 


where 


= 2 mone: 1/2 
oc ee ae ie 
an eevee a Nie + oy = nd, 
AX. — =. x 6 6 ed! end 
im te x x 
i ne eZ, 
te O 
X(£): complex frequency spectrum of the transmitted 
electrical signal. 
ee = cly, +nd ); speed of sound at transmit 
element with index n. 
bee See Gye mcd) eindgnittude Of Initial propagation 
vector at transmit element with index n. 
ow d complex weights of the transmit array. 
HYy(f,£ ;y): Random time-invariant space-variant 


transfer function of the ocean medium. 

NOt etatuis tS aderuncrlom Of beth indices 

Cendant, ier, 1S a function of the y 

eee acl ale weOumEme source Necatlon yo + nd 

iicmmaecel eur hoeatilon sy. +qd..-. This canY 

also be expressed by using the notation 

Hig, vo. ,0,q) expressing H, as a function of 
ne. © . 

wave number k_, direction cosine v_ and 

indices n and q. 2 


The above equations show how the medium transfer 


function is used in describing the acoustic field at a receiver 


Meecaelon as a fUnctEIOn Of Ehe acoustic field at a transmit 


Mmecation. 


The far-field directivity functions of transmit and 


i 


receive arrays are used to couple an electrical signal to 
the medium and vice-versa. 

The region of integration over direction cosine We 
is limited by. Sa This expresses the fact that the transfer 
function 1S not valid fer dinectton coo ie a yaies Ve approdemas 
ing 0 (close to a turning peéint). In most applicacivens en. 
limits of integration for both integrals can be further limited 


as ‘shown In Secei1on 


B. CONSTRAINTS ON THE MODEL 
In this section the inherent constraints on the model will 
be investigated. Also, some additional constraints will be made. 
1. Geometrical Constraints 
These constraints concern the transfer function which 
is only valid under the two following restrictions: 
a. The assumption made in the derivation of the transfer 


function, in order to validate a binomial expansion was [Ref. 2]: 


| {n7(y) - 1) /ve} | el (22s) 


Thus before using the computer model one should evalu- 
ate equation (2.23) for the chosen sound speed profile and 
array locations to make sure the inequality is obeyed. 

b. Absence of turning points. The eceunrenco et yceurnine 
points in a particular transmitter/receiver geometry depends 
on the initial direction of propagation, described by the 


angle By between the initial propagation vector and the 


20 


Y-axis. Using Snell's law we find that at the turning point 


depth Yq? 


Clq) = co/sin(8.) (224) 


The relation between the angle Ee and Nin 1S given by: 


< 
iI 


cos(6.) C2325) 


Combining equations (2.24) and (2.25) into an inequality 


which has to be obeyed to avoid a turning point results in: 


Cc 
Oe ee 1) eae (2.26) 
) 


Gy) 


and we conclude that obeying the inequality (2.23) also 
guarantees the absence of turning points. 
2. Sound Speed Profile 
The results presented in this thesis are based upon 
a sound speed profile which is assumed to have a constant 
gradient g and a random component of the index of refraction 


n, which is not a function of depth. Note that the equations 


R 
derived in Ziomek [Refs. 2,3] are more general and that the 
computer program described in this thesis can be easily modi- 
fied to handle a more complex sound speed profile. For purpose 
Beabcosuter Simulation, the random component of the index of 


refraction was assumed to be a zero mean, Gaussian random 


number with variance o* and was modelled by a random number 


Za 


generator. These restrictions make it possible to evaluate 
the transfer function with the closed form expressions given 
in equations (2.11), (2:16) and (2.19) "ter vane ieee et. 
ministic phase component and random phase component of the 
transfer function, respectively. In the calculation of the 
random phase terms of the transfer function, a different 
random number was drawn for each combination of transmit 
element depth io a and receive element depth acer 


thereby assuming that the random phase terms for the different 


transmission paths are uncorrelated. 


3. Frequency Spectrum Of the Treansmlmecd. fc cmiemecs 
Signal 


The transmitted electrical signal was represented by 
a finite Fourier series with fundamental period ie and maximum 


frequency tae HZ 


£ = KMAX (22-299 


oo 
max T 
O 
where KMAX denotes the total number of harmonics in the signal. 
This representation does not impose a severe restriction on 
the input signal, since every finite energy signal can be 
represented in the sense of a least mean-square error by a 
finite Fourier series. The expression for the frequency 


spectrum becomes: 


KMAX 
ee) ues ie C(£ eke ) ee er oes (2.23) 


Pa, Ou 


where fy is the fundamental frequency rs andwthice DC term 

is assumed to be zero. The complex Fourier series coefficients 
CL determine the amplitude and phase relations among the 
different frequency components. The integral in equation 
(2.22) now reduces to a summation and the expression for the 
real output electrical signal at a receiver element with 


indices 1 and q becomes: 


“y 5 . a 
y(t,i,g) = 2Ref (kf) 
k=l io nec 
q 
ae ae g Dict eit ae 
ne gn! OPI IK, MO) AZ) 
eel) 7/2 
fi 72 c_ exp [-jk u Ax, dv, du, exp{j2nk£t] }; (2.29) 
2 


O 


Zs 


TIi.), IMPLEMENTATION OF ieee epre 


The main task of the computer simulation model is to 
compute and file the electrical output signal at all elements 
of the receive planar array according to equation (2.29). The 
model was implemented on an IBM system/370 mainframe computer 
uSing FORTRAN. The program is computation intensive because 
the evaluation of equation (2.29) involves a double integral 
with nested summations. Therefore, much attention has been 
given to program efficiency. This section discusses the 
equations which were actually implemented and the overall logic 


of the program. 


A. NUMERICAL METHODS 

The model computes the real output electrical signal at 
all receiver elements as given by equation (2.29). Note that 
the double integral with respect to (wer st. saicecener 
cosines US and ve behaves like a frequency response or trans- 


fer function H(f,1,q), which depends om the position Gf iene 


receiver elements as indicated by the indices (1i.q). We define 
A £2 f f — 
Hien) = d (k V 1/q) 
of “Loa, ne (1) /2 nM no 


exp [-4k lexp[-jk_ (l-u2-v*) 7 2az) ree exp [-3k lacee 
-7 AY —7 —u - a7, Ee —jk AX. u u 
J*n%o qn : ane Vo! ion m I*n imo OO 


oe 
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Note that the factor L/e* in equation (2.29) was moved in 
front and replaced by the factor Wenn This 1s possible 
because the differences in sound speed at the different trans- 
mitter elements are very small for all practical cases and 
the effect on the amplitude of H(f,1,q) can be neglected. 
However, the dependence of the initial speed of sound on the 
index n of the transmit array will be taken into account in 
the evaluation of the phase terms as expressed by the different 
values of KO iP se guat rot (3.04). . 

Wetag the definition Of H{fi,i,q), the expression for the 
real output electrical signal at a receiver element with 


indices (i,q) becomes 


KMAX 
eae) = see Ret =} eae eco emit =) ) (3.2) 
k=1 : . 


ine computer program model first calculates H(f,1,q) 
as given in equation (3.1). Then, with the use of H(f,i,q), 
the output electrical signal is readily computed according 
memeacuation (3.2). 

The ocean medium transfer function Hu 1s possibly random 
and therefore difficult to integrate numerically so a change 
in the order of integration in equation (3.1) is made to 


achieve a more suitable form for implementation: 
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Fo a oe aus 
H(f£,1,q) = = dqd.H.(k_ ,v_.n,q) exo ljhs eee 
ce Seen pee ee ee 
Oo” 
1 (M1)/2 


c_exp[- jk, Uy AX. im OXP LO Ik (l- Us eye) Whee 


AZ)du_ dv (3:43) 
-l me-M-1)/2 ™ oo 


In order to improve efficiency, we concentrate first on 


the inner integral w.r.t. Uo. which we rewrite as 


i 


o ee 2 alee 
ae S,(£,u,)expi-k, [u 4X, + (1-ul-v")°/ ©Az] tau, (3.4) 
where: 
AR = S%> = 2S rd 
a s O x 
Note that 
~~ 
See Ue) Cc. exp i 7k area (375) 
X O =e ie eo Now. 


is the transmit far-field directivity function of a line array 
consisting of equally spaced point-sources as a function of 
direction cosine Uo and frequency f£. Now equation (3.3) reduces 


cor 


2 1 (wl) 


BCG i) a ) » Satan V.m,q) exp[-jk Vv AY] 
re n=- (N-1) / 4 


Q 


i 


= 2, 1/2 
f S,,(£,u,) expt jk fu Ax, fee us avs 


a AZ] }du, dv. (25) 
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In the case of simple amplitude weighting in the x-direction 
such as rectangular, triangular or Hamming weighting, we can 
@etarn closed form expressions for Sy (£,u,) - This improves 
program efficiency. For example, in the case of rectangular 


amplitude weighting 


Siac — un Ma 7 2 
(2g) = pare eae 
nee Dy x 
where Uy, is the direction cosine value to which the transmit 
far-field beam pattern is steered. 

Before presenting more details on how H(f,i,q) 1s calcu- 
lated, a brief overview on numerical and analytical integration 
techniques which were used is given. 

1. Integration Techniques 


a. Method of Stationary Phase 


This method approximates an integral of the form 


b 
ieee Ow mec (zexp [aig z) Idz (3.8) 
a 


where g(z) 1S a rapidly varying phase term compared to the 
slowly varying amplitude term S(z). Officer [{Ref. 10] shows 
Eiaeetaits Integral can be approximated by using the fact that 
the main contribution exists at the stationary point Z which 


is defined by 


The limits of the integral may be extended to infinity 
because only the region around the stationary point contributes 
to the integral. Assume that around the stationary point, z 


is given by 
7 i ee, + 2 (3.108 


The approximation of the integral is then given by 


z= 5 : TI 
t = S(z.) 2a ala (2,,,) lexplilg(z..)4q]] (3.299 


under the condition: 


Logan (2. ) | 
scl ee C3 Lea 
cutee 


where the sign in the phase term of equation (3.11) equals 
the sign of the second derivative gee 

b. Numerical Integration Using Newton-Cote's Formulas 

This standard technique [Ref. 11] is the method 

used by the program to perform a numerical integration over 
a subinterval of the total region of integration. Subinter- 
vals can be defined by dividing the total region of integra- 
tion into equal parts or by a scheme as described in the next 


section. Each subinterval is integrated by evaluating the 


Simpson's rule estimate: 


(s + 2s, + 485. + 28. + 45, + ... + 28 + s_) 
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where Sy a S denote values of the integrand evaluated at 
different points of the subinterval, and the distance between 
these points is given by the step size h. The original step 
size h is equal to half the size of the subinterval to be 
integrated. The step size is then halved to obtain a more 
accurate estimate which takes more points into account. This 
is repeated until a tolerance criterion is met. Using succes- 
Sive Simpson's rule estimates, one can do a Richardson's 
extrapolation [Ref. 12]. This extrapolation uses the last 

two estimates by Simpson's rule, §S 


D and Shr where S., uses 


2 


half the step size of § Both estimates have a global error 


1° 
©f Order h*, The extrapolated value becomes: 

Extrapolated value = S.5 + (S. ~S,)/195 (S74 ) 
with an error of order Bice If two successive extrapolated 


values are within tolerance, the last extrapolated value is 
taken as the value of the integral for that subinterval. This 
scheme creates a semi-adaptive integration procedure because 
different subintervals can require different step sizes. It 
could be further enhanced but would then require more overhead. 
Another possibility would be the use of a Gaussian quadrature 
rule which needs less function evaluations for the same order 
of error but has the disadvantage that the points are not 
identically spaced aw the scheme cannot be made adaptive in 


a computation efficient manner. 


Zo 


c. Numerical Integration Using the Euler Summation 

A problem in evaluating an integral like equation 
(3.4) is convergence. The real and imaginary parts of the 
integrand of equation (3.4) are oscillatory. Successive 
positive and negative contributions tend to cancel, but each 
contribution itself is not negligible and the integral behaves 
oscillatory as a function of its limits of integration. There- 
fore, one has to integrate over a relatively large region to 
obtain convergence. This is in conflict with the desired 
computational efficiency. 

Squire (Ref. 13: pp. 172-173) describes vane mes 
gration procedure for oscillatory integrands which splits the 
range of integration into parts using the zero crossings Of 
the integrand as points of separation. The individual subin- 
tervals can be evaluated with the.use of standard methods 
(as described in III.A.1.b) and summed to form the value of 
the integral. A technique like the Euler summation [Ref. 13: 
pp. 172-173] can be used to improve convergence of this sum- 
mation. The Euler summation obtains from an original sequence 


Aa al See el 
ars 2a ent eT 


a second sequence 


and a third sequence 
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Cc) = b. + by. Cy = b, tbo, 7, Cf = boy thy 
and so on. It can be shown that the sum of 
dl ui I al 
7 4] + 7 Dy earl ee, or 
is the same as the sum of the original sequence. Convergence 


however is much improved which cuts down on the number of 
subintervals to be integrated and improves efficiency. 

This method works well on integrals where successive 
contributions from subintervals to the integral are of decreas- 
ing magnitude and of alternating sign as illustrated in Figure 
2. To use this method it must be possible to solve for the 
zero crossings of the integrand, which generally requires a 
closed form expression to solve for the roots of the integrand. 

Pee OD lication Of Integration Methods 
The methods introduced in the previous section are 
used to define a procedure to calculate H(f,i,q). First, 
the integration w.r.t. direction cosine Uo is considered, 
Widen 1s given by 
: 2 2, 1/2 
O 


if Pea expi- gk lu Axed lou —v AZ] }du, (3.4) 


Both the method of stationary phase and the numerical 
method (Sections TIZ.A.1.b and IJII.A.1.c) using the Euler 


Summation are implemented. To apply the method of stationary 


Sill 


INTEGRAND 


Gt Z ) 
end SUBINTERVAL 
Ist SUBINTERVAL 
Bato ee Oscillatory Integrand Showing Decreasing 


Contributions to the Integral 


BZ 


phase we define in accordance with equation (3.8) 


ae ely 2 
g(u) = kK, l4x;u + (1 US ae AZ] (35) 


The amplitude term S,(f,u,) is the directivity function which 
is indeed slowly varying w.r.t. g(u_). We find the sta- 


Memary point UOSP by setting g'(UOSP) equal to zero: 





; =. Se oe - 
g'(UOSP) = k 4%, +k 421 US a eg = OCs .ccl6)) 
leading to 
Z 
AX. i ae es 2 
UOSP = Az ! ) (324) 


ene ey Az) 


For the values of the second and third derivatives one finds 


i 7 me 8/2 
g (uo) = k 4201 us Le ‘el ad es 
Wwe ae an Z - 2 ay a2 _ Z. 
ua (ey eer W(t ny 2) aw) (3.19) 


We now form the ratio of the third and second derivatives 


at Uo equal to the stationary point UOSP: 


Z Jae MN ee 
Gm OOr = _ 3AX eee i Seon 
a CUCsSE) 4h O Age i 


The condition in equation (3.12) becomes 
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Z 


AX. il + (AX, /AZ) 
c (— =>) ——— ace (22 21) 
AZ Z 
dk aN. 


from which we conclude that the above condition is obeyed 
only when AZ >> AX., which means that the method is only 
valid for relatively small cross ranges XTX Since 
AX, <= it ae a ie 
a 16 O x 
The program implements the following procedure to 
compute the magnitude and phase of the integral 772. us 


[see equation (3.4)], denoted by IWRTUO. 


1) Calculate the stationary point UOSP 





AX. Jt -ve 
UOSP = AZ oS (3 is 
1 + (AX, /A2Z) 


2) Calculate the value of the second derivative DRV2ND 
at the stationary point: 


DRV2ND = kAZ(1 -UOSP* - v*)~°*/*(1 - ve) (3. 2h 
3) Calculate magnitude and phase Of “Ehe -intecwal: 
MAG = S(f£,U0SP) /21/DRV2ND (3.23) 
PHASE = -k  [AX, UOSP + (1 - vosP* -ve) re +7 


Note that the integral is a function of ee and the indices i 


and n through KO and AX, . The dependence of the magnitude 
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Smee index nis negligible, thus DRV2ND is calculated for 
Ewer USing Ko. We will denote the integral by IWRTUO(v.,k_,1). 

The numerical method using the Euler summation can 
also be applied to the integration w.r.t. Uae This numerical 
method is implemented to verify the result of the stationary 
phase method and extend application of the program to larger 
cross ranges. The applicability of the method can be shown 
by separating the integrand of equation (3.4) into a real 


and imaginary part 


al 


2 Zeya 
a (S(£,u,)cost-k, [u Ax, + (1 UE ae Le AZ] } 


(nee 5) 


2 


stu.) simi-k lu ax, + (i mu? -v2) 1/2 


AZ] }) du, 


As shown above, the point of main contribution to the 
moeegral in equation (3.25) is the stationary point UOSP. 
Away from this point the integrand becomes more and more 
oscillatory and consequently more difficult to integrate 
numerically. One would therefore like to use a scheme in 
which the size of the subintervals is dependent on the behavior 
Of the integrand. The method described in III.B.l.c uses 
the zero crossings of the integrand as points of separation, 
thus making sure that the integrand within a subinterval is 
"well-behaved" and can be efficiently integrated by standard 
numerical methods. This division into "well-behaved" subin- 
tervals could also be obtained by using the relative extremes 


aS points of separation. The integrand of equation (3.25) 
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has closed form solutions to the zeromeressings anger e tae. 
extremes of both the real and imaginary parts. The zero 
crossings of the imaginary part of the integrand are found 


by solving: 


=e JNK ez Can nt ; n intéger (3.26 
n iL oO 1g) O O 
or 
~AX= (am) AZUL 2 ee ea eee 
ale O PP nh 
CR eS (3 oe 
Kien eee Zee) 
i file 


These also are the locations of the relative extremes 
of the cosine term appearing in the real part of the integrand. 
In order to avoid integrating the real and imaginary parts 
of the integral separately, the same subintervals are used 
for both the real and imaginary parts of the integral. Points 
of separation for the subintervals are the zero crossings of 
the imaginary part of the integrand, which coincide with the 
relative extremes for the real part of enew im cea mane 
starting point for the integration, define the first subinter- 
val using the point of main contribution, i.e., the stationary 
peoine UOSr= 

Fig. 3 shows a typical behavior of the real and imagin- 
ary parts of the anieege ana around the stationary point.” Tine 
figure also indicates the position of the subintervals. 
Integration of the successive subintervals (from UOSP outward), 


will form alternating terms of decreasing magnitude for both 
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"MAGINARY 
== - REAL 


Ud 


NZXIT 


Eaves 


Meg] 5 Typical Behavior of Integrand in 
Peon ec netmewctattonary Pont 


OF) 


the real and imaginary parts of the integral, making it 


possible to effectively speed up the cConvergencams meinem 


complex Euler summation. The decreasing magnitude is caused 


by the non-linear behavior of the phase term and the decreas- 


ing amplitude of the integrand which is given by the transmit 


far-field directivity function S,(£,u,)- The program imple- 


ments the following procedure: 


1) 


2) 


=) 


Calculate the value of the stationary point UOSP 
according to equation (3.17). This value of u_ also 
denotes the maximum (negative) value of the phase 
coun: 


Calculate NZ the maximum integer multiples of mT 
contained in the maximum (negative) value of the 
phase term (see Fig. 3). 


Be [UOSPAX , +(1-U0SP*-v*) */*AZ] 


NZ = INTEGER ( | —+~——-_=____~___°______}) (3.289 


Note that NZ denotes the zero crossings of the imaginary 
part of the integrand nearest to the stationary point 
UGSE. 


Calculate the values of u, £or the Zeromcervoss ings sce 
the imaginary part of the integrand denoted by NZ 
Using equation: (3520). 


=AX, (NZ) £AZ [ Cle ee = Ne ae aa 
UOZ, 5 = = O 2 A 
k n AX, +42 ) 
(3 225) 


Integrate the subinterval between these two values 
OF UOZ. 


Decrement NZ and calculate the values of UOZ), 21 
denoting the next zero crossings of the imaginary 
part of the integrand, by using the formula given 
in 3) above. 
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5) Define two subintervals, one on each side of the 
Stationary point, from zero crossing to zero crossing 
for the imaginary part and from extreme to extreme 
for the real part of the integrand. 

6) Integrate both subintervals and add the contribution 
to the total value of the integral using a complex 
Euler summation. 

7) Check for convergence by comparing the difference 
between the last two values of the Euler summation 
with a tolerance. If not convergent, go back to 4) 
and repeat the procedure. 

One more note can be made concerning the numerically 
@alculated value of the integral w.r.t. ua Poe Cle terent Ko: 
Recalculation of the whole integral for different values of 
n would be very time consuming. We therefore introduce a 
correction term to obtain values of the integral for n #0, 
based on the value of integral for n= 0. This correction 
term is only of importance to the phase of the integral, the 
effect of different values for KO on the magnitude is 
neglected. For n = 0 the stationary phase method results in 
the following expression for the phase: 


2 


PHASE (n = 0) = -k, [AX,UOSP+(1-UO0SP*-vg)°“/*Az]+ Z (3.30) 


a 
Brekhovskikh [Ref. 14] demonstrates that taking into 
account higher-order terms in the method of stationary phase 


does change the amplitude but not the phase of the resulting 


approximation. This shows that the point of main contribu- 
tion, i.e., the stationary point determines the phase of the 
integral. If the method of stationary phase is no longer 
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valid we can write a general form for the phase based on 
equation (32207r: 
PHASE (n=0) = -k [AX , UOSP+(1-U0sP*-v*) 7721 + e oe 


(See 


where Ad represents some unknown correction term. This relies 
on the fact that the stationary point UOSP is always very 
close to the point of main contribution. Based on this, the 
phase term for n # 0 can be written as: 


2 


ae = _ 2 | Z ae 
PHASE(n) = [k + (k kK) 1 (4x, UoSP+ (1 UOSP vo) AZ)+7 + do (3.32) 


where 


(k_-k_) [AX,UOSP + (1-uosp?-v2)*/? 
n “o it . 


AZ] (333) 
represents the correction term which is readily computed for 
different values of n. This phase correction term has to 

be added to the phase for n = 0 to obtain phase values for 
n#0. Note that this correction term is only important for 
larger ranges AZ, assuming that AZ >> AX. This can be seen 

by calculating a typical value for Kee where the element 
Spacing is taken as \/2 and the sound speed profile has a 


constant gradient g: 


Q|A 
ka 


oT al 58 = Me = = ae) SS nN =) BG WOM eens (3.39 
ro E GC 


1g O 


O 
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where we used Chae 1500 m/sec and g = 0.017 1/sec. An upper 
mae fOr the correction term is (K-K9) AZ ate ere 

Having established two methods to implement the calcu- 
Beerom Of the innér integral w.r.t. direction cosine Uo: the 
Sleer integral w.r.t. direction cosine Ve 1s now considered. 


The integral can be written as [see equation (3.6)]: 


il N= dij 2 
Gumus. 7ve,n,oexol=7k.v_AY 
a PesN=1) 72 eee lO io 


J IWRTUO(v_,k pe CA 


qn n O 


W335) 


The ocean medium transfer function values are calculated with 
the use of equations (2.11), (2.16) and (2.18). Note that the 
deterministic and random phase terms pertain to the case of 

a sound speed profile with a constant gradient. In applying 
mmese formulae, we have to take into account that the source 
gepth for a transmit element is not ee ut y,tndy - This results 


in the following formulae which were actually implemented: 


An = (eae) mee Base 
n a n 
6 (Gacy) = »=—[—(——-; - l) + AY__] (3234) 
MD 2V, g BG Sie) qn 
sa Ch oat 
Ou ined) a vo' 9 ° TNR on Sty ads) (Seas 6) 


We note that both the amplitude and phase terms 


depend on the value of the wavenumber KO and direction cosine 
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a The expressions within the square brackets for oMD and 


Q are only dependent on the geometry and can be precomputed. 


MR 
A random number generator was used to generate a series of 


N(0,1) distributed numbers to be used for n in equation 


NR 
C323 ae 

Because of the complexity of the total phase term in 
the integrand of equation (3.35), a closed form approximation 
could not be found. The numerical method used is also based 
on the procedure given in III.A.l.c. The integral is divided 
into subintegrals in such a way that the subintegrals are 
all "well behaved." The problem is to find a general Foun 
which approaches the phase of the integrand. This form can 
then be used to calculate the points of separation in order 
to define the subintervals. We will proceed by separating 
the phase term of the integrand into a major part and a minor 
slowly varying part. Combining the phase term of the integral 
Wo ce un given by equation (3.32) with equation (3.35) 


results in the following expression for the integration w.r.t. 


Vos 
O 
1 (N-1)/2 
as { ee dB Vg a expt—jk Vv, (ad — nd.) +7 +A} 
q 


a q : 2 ly ; 
exp{ j(k, iy [AYv +X , UOSP+ (1 UOSP Vv.) AZ} } 
. Die? 
| IWRTUO(v_,k_ 4) jexp{—jk_ [v A¥+U0SPAX, + (1-UOSP v5) 44] }dv. (3.39) 
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where AY = iss Ree If we let 


(N-1) /2 
M(v_) = y dH (k_,v_,n,q)exp{-jk_v_(qd' -nd_) +++ Ad}- 
O n==(N=1) /2 nM He @ ne @: ay y 4 
exp{-j (k_-k_) [Ayv +X. UOSP+(1-UosP*-v2) “7az] } + | IWRTUO (v ewe Dy) ques 5 ie) 
A, Oo Ovvicr O Ou i ° 
and 
a(v_) = -k_[AYv ec pa eCIsUCSP——7-) “AZ 
O O O 1 O 
(2220) 
a(v_.) = -k_[AYv eee 7 Vix ae 
O O O O 1 
Then the integral given by equation (3.39) becomes: 
1 
; f X(vo)expl[ja(v)] dv, 
g 
where X(v_) is slowly varying compared to exp{ja(v.)}. The 


phase term atv) can now be used to define subintegrals which 


are "well behaved." This is accomplished by taking the roots 
OL sin(la(v.)) eis cos(a(v.)) as the points of separation between 
the different subintervals. Contributions from the subintegrals 


are in general of decreasing magnitude and alternating sign, 
so the Euler summation can be used to effectively een up 
convergence. 

The implemented procedure uses the roots of sin(a(v)) 


as the points of separation and is similar to the one already 


43 


described for the snereqration w. race Ua: Tie pointe of Gan 
contribution is approximately where the phase term av) is 
stationary. This point is used to begin the numerical inte- 


gration. The stationary point VYOsr as givens, 


VOSP = AY¥/(RO, + a¥*)7/* (3.43) 
which is a solution of a'(VOSP) = 0 where Ros is the radial 

distance from the center of the transmit array to a receiver 
element with index i and is given by 


Hee 
R.. = (AX. ietean) (3.44) 


The points of separation for the integral can be found by 


DOR nee Eee 
—AY (nm) Ries ay Ot Re ieee 
_ eo ia 
; k (RO. + AY”) 
Oana 
Another method using a brute force technique was used 
to validate the above described procedure. This method 


integrates w.r.t. direction cosine Vee between user specified 
limits of integration. It divides the range of integration 
into subintervals of equal size which are integrated by a 
standard routine as described above in III.A.1.b. Both 
methods agreed numerically although the method using the 
unequally spaced subintervals and the Euler summation was 


much more efficient in terms of computation time. 
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Be. 


PROGRAM OUTLINE 


The flowchart in Fig. 4 gives an overview of the program 


ieegre. fhis logic is controlled by flags which are read 


from the input data file together with all program parameters. 


The 


program parameters determine: 


ReaeiacOMmo ety. <0) yo eee, Voy, 2, M; N, M*, N', 
- ee ce O O O fe te ir 

x i x y 

Speed of sound profiles: Cor Gr O- 


Beam steering angles for the transmit array and the 
Frequency to be used for the beam steering calculations. 


ie cundamental period T, of the electrical input signal 
and the total number of samples L in the electrical 
output signal. The program determines the sampling 
poatocd 12 — ie7] and calculates the output electrical 
Signal at time instants 


—(L-1)T /2, Raia, | Cue. (L-1)T /2 


The values of the frequency components in the electrical 
input signal and their complex Fourier series coefficients 


Cy 

Tolerance values ERRORU and ERRORV for the numerical 
iteeqrations W.rst. u. and ve The user should validate 
his choice for the tolerance°values by Comparing them 


against the absolute values calculated by the program 
for the integrations. 


Processing is done in accordance with flag settings and 


pielows for a number of options: 


The frequency response H(f,i,q) can be calculated or 
read from a previously created file. 


The calculated frequency response H(f,1i,q) can be 
filed for later use and/or printed. 


The output electrical signal can be calculated and 
filed for further processing by other programs. 
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Read input 
from file 


= 
= 
SE 






Calculate 
Mitts Ls a) 









Read 
MCT, lL, cl? 
from file 









Write 
oo CL th 
to file 












Calculate 
output signal 
Y(t) 


EE 







Write 
¥ Ceri) 
to file 


Fig. 4 Flowchart of Main Program Module 
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Most important numerical calculations required by the 
model are done within the module CALCHF. The module CALCHF 
calculates the frequency response matrix H(f,1i,q), i.e., the 
frequency response for all receiver elements (1,q). This 
weeny involves evaluation of the double integral w.r.t. u 
and Vo: It has two versions dependent on the integration 


technique used to evaluate the integral w.r.t. uo? 


CALCHF 1 Use stationary phase method to integrate 
yee er Bl 


CALCHF 2 Use numerical technique to integrate w.r.t. uo: 


To avoid recalculation and improve program efficiency 
maewouter integration w.r.t. Vs 1s done simultaneously for 
all elements with the same index i. This is particularly 
important when using the more time consuming numerical method 
to integrate w.r.t. Uo: This can be done by separating the 
Summation inside the integrand into a part dependent on gq 
and a part dependent on i. Recall that equation (3.35) 


gives the form of the integrand as follows: 


Mi) /2 
a on ;. 4 Hy (kV 1@) expt Jk VAY JIWRIUO (Vv, 1K, i) (3.46) 


SO given a value for IWRTUO for a given index i we can evalu- 
ate the integrand for all values of g, thus avoiding recalcu- 
lating the most computation intensive part. This scheme does 


result in more complex coding. The program has to do N' 
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integrations and keep track of N' Euler summations simul- 
taneously. It does, however, also reduce some overhead. 
Convergence checking on the integration w.r.t. ee is done 
only on the two outermost elements with indices gq = (N'-1)/2 
and gq = <(N'=-1)/72. 

The calculation of those terms that multiply IWRTUO(K Vo ,1) 
in the summation given by equation (3.46) is simplified by 
precomputing part of the phase components of the ocean medium 
transfer function given by equations (3.36) and (3.37). The 
program does this in the subroutine HMWKB and stores the 
resulting phase integral values in an array. To obtain a 
particular value for the phase, the program only needs to 
multiply wath theweaceor eal Especially in the case of a 
more complicated sound speed profile, the computation of 
the phase integrals can be involved and precomputation im- 
proves efficiency. The stated time-invariance of the model 
also forces precomputation since the random numbers used to 
model the phase term OMR can be drawn only once. The random 
numbers were obtained from the IMSL routine GNNML which was 
used to generate a series of N xN' random numbers with distri- 
bution N(0,1). From these random numbers the different random 
phase terms Our for all combinations (n,q) “aneruecadily con 
puted using equation (3.37). Figure 5 shows a simple flow- 
chart of the CALCHF module. 

The innermost block which actually calculates H(f,1,q) 
1S shown in more detail in Fig. 6. The diagram shows how the 


subintervals are defined by the program. The procedure Toltevee 
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CALCHF 


Call HMWKRE 


Calculate 
phase integral 













DO 
for all frequency 
comoonents 


DO 
for allen 






Calculate 
ruGhe tee 
ron sal a 


e 
S 


Pelee plewecwomemon Module CALCHE 
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Calculate 
V@SP, NZ 


Set: 
V@BEG1 
V@BEG2 


Using NZ calculate 
zero crossings 


V@END1, V@END2 
(Eq. 3.45) 


Define intervals: 
right: V@BEG1 — V@END1 
lett: V@ENDe ~ V@BEGe 


Integrate 
subintervals NZ = NZ + 1 
VOBEG1 V@END1 


V@GEGe VJEND2 








Add contribution 
to Euler summation 





convergence 


VOSP: Approximate point of main contribution of the integral. 

NZ : Maximum negative integer denoting the multiples of 
contained in the maximum negative value of the phase. 

Note: The “Integrate Subintervals” block makes a Subroutine 
call to evaluate the integrand given by Eq. (35.46). 


Fig. 6 Calculation of Subintervals by CALCHF 
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memcaescribed in IIIT.A.2 for the integration w.r.t. ie and is 
Similar to the procedure used for the numerical integration 
ol an US iM Sect 1Ome inl Av. 

Implementation of the integration over a subinterval using 
the Newton Cote's formulas is straightforward and will not 
be examined in detail. The framework for this part of the 
feeegGam Can be found in Squire [Ref. 13: p. 281]. 


The Euler summation is used to speed up convergence of 


the numerical integration. The objective is to form the sum 
ee, 
2 4 8 16 

fem the Original sequence Ay 1A51 Az Ay, very AL which are the 


contributions from the subintervals 1 through n, where 


Cc =o ae ie 


an oj k . aa? 


k k k+1’ 
The algorithm used is described below and uses an array 


SOBF to store the coefficients ane b c dq etc. 


M=1" ne-2° “n-3’ 
The size of this array determines how many contributions can 


be handled. 


DIMENSION COEF (2,SIZE) 
REAL EULSUM, FACTOR 
INTEGER NEW, OLD 

ce AGE coneri bution N 


COEF (NEW,1) = N_ th CONTRIBUTION 


ay 


DO 100" J 22a 

100 COEF (NEW,J+1) = COEF(NEW,) J). =3¢@Gri. ae 
FACTOR = PAGIOR 32 
EULSUM = EULSUM + COEF (NEW,N) /FACTOR 

E Prepare fOr Neste oir outro 


SWAP( NEW, OLD ) 


The program implements this algoritm in complex form to 
perform N' summations simultaneously. 

More details on the actual programming of the discussed 
numerical methods and flowcharts, data structures used, etc. 
is considered applied programming and will not be discussed 


in this thesis. 
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Pe COME Emo LATED DATA 


This section describes the validation of the model. It 
presents and interprets the data generated by the computer 


meaqel for a number of test cases. 


Pwr ROE NTATION OF COMPUTER SIMULATED DATA 

Data generated by the computer model consists of the 
frequency response matrix H(f,1,q) and the output electrical 
Signal at each element y(t,i,gq). Interpretation of the data 
is done directly by examining the frequency response H(f,i,q) 
and by the use of a three-dimensional DFT program to calculate 
the Fourier transforms of the computer simulated output elec- 
trical signals w.r.t. time and space. The program is based 
upon 3-dimensional DFT beamforming for planar arrays as des- 
cribed by Ziomek [Ref. 4]. First a brief overview of the 
Smee, trom the 3-D DFT program iS given. 

The program calculates the 3-dimensional DFT of the 
electrical output signal from all receiver elements and is 


given by 


M! Nit le 
ec,r,S) = } ) )  ¢ diy (&£,m,n) exp [-j2mg2/L] 
m-M' n=-N' 2=-L' 


exp [j27xm/ (M+ Mz) Jexp[j27sn/ (N+NZ) ] (a |) 


DS 


Wikeiee:: 
L: total number of time samples in one fundamental 
period aA of the signal. 
L' = (L-1)/2 
T.: sampling penued 


M: total odd number of elementsS in the x=direce ten 
in the receive array 


MY = (Mei 2 


N: total odd number of elements in the y-direction 
in the receive array 


N' = (N-1)/2 
da. ee Spacing in x- and y-directions, respectively. 
Ca qd: complex separable weights for the array 
given by 
cS a jexptj6} 
ci = b_ exp{j$,} 


q,r,S: binnumbers of the DFT which are related to 
values of frequency £, direction cosine u 
and direction cosine v, respectively. 


y(2,m,n): output electrical signal at time instant 
ATL at element (m,n) of the receive array. 


MZ, NZ: Number of zeros used for padding to increase 
the resolution of the direction cosines u 
and v, respectively. 
Note that the symbols used to characterize this planar receive 
array are not the same symbols used to describe the receive 
array in the computer model (i11:-A\l). Réectangulan ampli cue- 


weighting without beam steering is used for all the results 


presented in this chapter. 
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In II.B.3 we made the constraint to represent the transmit 
electrical signal by a finite Fourier series with maximum 
frequency KMAX*f.. Thus, to avoid aliasing, we must sample the 
received signal with a sample rate equal to or greater than 
the Nyquist rate: 

t. a 2 KMAX£. ; To a 2 KMAXT . Car 2) 


The total number of samples L to be taken in order to avoid 


aliasing must obey the following inequality: 
ioe ee eee (4.3) 
The sampling period is determined by: 


To = TO/t (4.4) 


The program evaluates the 3-dimensional summation in 
three steps. The first step is to perform a DFT w.r.t. time 
at each receiver element (m,n). If the sound field incident 
upon the array is a general plane wave, then the received 


electrical signal can be expressed by 


{| 
WQ 
a 
| 

l 
a 


7 (ely 1a) 


where t_ = (u_md. + v_nd_)/c denotes the time delay between 
O @ 2s oy 
the signals at the different elements. The DFT w.r.t. time 


becomes 


5 


i Fas 
= 2 q* 
Yo(q,m,n) Cia — g (kT. to) Wy (45) 


where we” is the factor exp{-j21mq2/L}. Using the above stated 
fact that the received signal is given by a finite Fourier 
series, we obtain 


Yo(q,m,n) = Con beg exPi = demgt (uma yond 7s (4 Gy 


A normalized expression is defined by 


You (4-m,n) a Yo(q,m,n,)/(c db) 
= cexpt-j2ngt (uymd, +vind)/e} (ae 
So You (G/m,N) allows us to estimate the amplitude of the fre- 


quency components in the signal at each element (m,n). 


The next DFT calculation is: 


mee 


Yo(q,xr,n) = ) Yo(q,m,n)exp{j2mmr/(m+mZ) } (4.8) 
m=-M' 


where we have padded with NZ zeros to increase resolution. 


Substitution equation (4.6) into equation (4.8) leads to 


Yo(q,r,n) = Led, exp{-j2n qi,v nd /e} 
M! ooh 
aie, cjexpt-j2ngf umd /chwWi vs (4.9) 
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and a normalized expression can be written as 


Yon (4,4) = ae exp{-j 21g See S/S 
M ’ M! 
ec expt en quence /ciWo 7 } a (4103 
et On 20 x M+MZ on i 
Mol m=-M 


mgs EXPressioOmM 1s directly proportional to the Fourier 
series coefficients a and to the normalized far-field beam 


pattern at a frequency q° ft, Omar liner dcrray Comsisting Of M 


elements which is lying along the X-axis. The planar array 
has N such line arrays as indicated by the index n. If no 
beam steering is done (Cc. = ant then Yon (Geb -n) has a maximum 


approaching a Pewee won Tuner related tothe direction 
cosine value which is closest to Ua: The DFT bin number r 


meecelated to the direction cosine value by 


uo = rA/[(MHMZ) a, ] ; h = c/qf, (4.11) 


This value of u is an estimate of Uo - Therefore, for each 
frequency component of a plane wave propagating in the wee 
direction, the expression given by equation (4.10) can be used 
to estimate the magnitude of the frequency component and the 
direction cosine Uo: An estimate of v.~can be obtained in the 


O 


same way by evaluating 


on 


Ce oe oF exp{-j27 qf, ujmd, /c} 


Ne ie 
a ,7on 
) . diexpi a 2mme Evo nd /clw Vy, La pL (aa 


een N+NZ/ 

The above mentioned estimates of Cat US and Vv, are provided 
by the program in the form of printouts of You (Gemn), 
You (G,45/n) and You (GeM,S) . The printouts are done for user 
specified values of m and n. No padding with zeros is done 
for the transform w.r.t. time. Padding for the spatial trans- 
forms is done to obtain a user specified step size in direction 
cosine u or v. Graphical output for estimation of un and va 
was also obtained. The complete three-dimensional transform 
Yo(q,r,s) was used to produce a three-dimensional plot to aid 
in interpretation and for illustrative purposes. The esti- 
mates of Uo and Vo can also be used to calculate estimates 
of the spherical angles denoting the direction of the sound 
Source as measured from the center of the receive array. 


These angles are also referred to as target bearing and eleva- 


tion angles. The estimates are given by 
py. = tan l(v fu) (4. Weep 
O oro 
5 _ alee 2 ee 
6 = sin (tue + ve ) (4.14) 


nN 


where Ws and ee are estimates of the target bearing and ele- 


vation angles, respectively. Note that an error in the 


58 


Seeeilactol Or Cither u, Or Vv Or both affects both the target 


O 


bearing angle and target elevation angle estimates. 


Eee LEST CASES 

Most test cases run by themodel are kept simple to aid in 
interpretation. All cases are run first in a deterministic 
homogeneous medium so that results could be easily interpreted 
and compared to known methods and approximations. All test 
cases are checked to obey the criteria reguired by the trans- 
feos function. 

The values of the basic parameters such as source depth 
Me? gradient g and reference speed of sound Cc, were chosen in 
order to model the wave propagation between a source located 
on the SOFAR channel axis and a deep receiver. At moderate 
latitudes from 60°S to 60°N the speed of sound on the SOFAR 
axis ranges from 1450-1485 m/sec in the Pacific and from 
mo-1500 m/sec in the Atlantic [Ref. 15]. The depth of the 
SOFAR axis is typically between 1000 and 1200 m at moderate 
latitudes [Ref. 15]. The gradient below the SOFAR axis is 
given by g = 0.017 1l/sec [Ref. 19]. Values for the standard 
deviation o of the random component of the index of refraction 
mange from 0.25 x1074 [Repomemi sees lito: 1.0 Sige: [Ret 16:1). 
These values were used in choosing representative parameters 
for the following test cases: 

1. Homogeneous Cases 

As a baseline test case we define the deterministic 


homogeneous case HMGI1: 


De 


Co 1475 m/sec. 

g = 0: constant speed of sound. 

o = 0: deterministic medium. 

M= 11, N= 11: all by 11 element transmit array with 
rectangular amplitude weighting in both the x 
and y directions. 

M' 59, N' = 5: a5 by 5 element receive array. 

od OF m K = o0e0 m 

O e 
we LEO OO rs 0 om ss ZOCOR Or = in 
Z OF m Zo = Li s2 20 on in 
O ie 
= qd'=q' = ° } 

d.. a Ce oe O.7379 Ms: Spacing eequal to Amin’ 2° 

0 SO, ee = 90°: spherical angles denoting the line 
of sight between the center of the transmit array 
and the center of the receive array. 

u 0, Ve 0.5: line of sight direction cosine values 
from the center of the transmit array to the center 
of the receive array. 

Uy, = 0, Vb - 0.5: direction cosine values used to steer 
the main lobe of the transmit far-field beam pattern 
towards the center of the receive array. 

KMAX = Single frequency component. 

E LOQOFRZ pera 10-001 cece 

O O 

Cy a, = 1.0: complex Fourier series coefficient equals l. 

Ll = 3: 3 times samples per fundamental period TO: 

Based on the definition of test case HMG1l we define 


the following 


HMG2 : 


Cases. 


Same as HMG1 except 
= 200 mm, Z = 1720-4648 
1 iG 
os uy, = Ore i i — ee 
= 30.66°, w.= 78.69°. 
ie ia 
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HMG3: Same as HMGl1 except 
ee = lO Omri, Cie soo .0246) m 
HMG4: Same as HMG1 except 
eo 25007 m; 2a 2598.074 m 
HMG5: Same as HMG1 except 
= Oe = ee do =O oo, Sm 
x y 


x BY 


fees 2000 SZ, i O05 sec. 
O O 


HMG6: Same as HMG1 except 


m=) 226.7425 ™ Zee 224 eS om 
ig ie 
o> Ve Omroreaia ea O25 
= ° = ° 
6. 52324; v B92 Bie 


HMG/: Same as HMG1 except 
(eee =a = d' = 0.184375 m 
Xx y x y 
£ = 000 Hz, ae = 0.00025 sec. 


HdMG8: Same as HMG2 except 


KMAX = 5: five frequency components. 

fee = AUC) Oat Z, ., ea 700 lesec ; 

O O 

1 =5000nH 2s 

max 

fh = 3000 Hz: frequency used for steering the 


transmit beam pattern. 


Gee — cee — Gd — d' = W47S/72 £ = 0.1475 m 
sx y x y max 
ome ae — | Os; kee 1 2,72..,5: all complex Fourier 
k k : ‘one 
series coefficients equal l. 
L= ll ll time samples per fundamental period To: 
2. Deterministic Inhomogeneous Cases — 


The following deterministic inhomogeneous cases 
incorporate a sound speed profile with a constant gradient. 


The definitions are based on those used for the homogeneous cases. 


be 


INHMG1, INHMG2, and INHMG5: Same as HMG1, HMG2 and HMG5 except 
g = 0.01721 /sec: 
INHMG8: Same as HMG8 except 


g = 0.017 1l/sec. 


3. Random Inhomogeneous Cases 
The following cases with an index of refraction com- 
posed of a deterministic and a random component are defined: 
RINHMG1, RINHMG2: Same as INHMG1, INHMG2 except 
6 = 10 alone 
RINHMG8: Same as INHMG8 except 


aloe eae 


C. RESULES 
This section will present and discuss the results from 
the computer model runs on the test cases as defined in Section 
IV.B. Some cases were run using both methods of integration 
Wie ey aie. uo: The two methods showed no significant differences. 
First, some trivial checks were made uSing the various 
homogeneous cases. Note, that all test cases have the same 
line of sight direction cosine v.- For the homogeneous cases 
with a single frequency component of 1000 Hz the amplitude 
and phase terms of the ocean transfer function are equal and 
given by: 


= ~l/f2 Z 
oe : 8 = 0 (4.15) 
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Since those cases have the same ocean transfer function 
values, the difference in the magnitude of the frequency 
response H(f,1i,q) for the various cases should only depend 


on the total range given by 





Oi aeE (4.16) 
and the transmit far-field beam pattern. Cases HMG1, HMG3 

and HMG4 also have the same line of sight direction cosine 
value UL, SO the influence of the transmit far-field beam 
pattern is the same for these cases and magnitude differences 
between the frequency responses can only depend on the total 
range and should behave in accordance with the expected 
spherical spreading for a homogeneous medium. Table IV.1 

shows cases HMG1, HMG3 and HMG4 together with their total 

range and the magnitude of the frequency response at the center 
element of the receive array, i.e., |H(f,0,0)|. We conclude 
that the magnitude indeed falls off as 1/r, only the total 
range is of importance. 

Table IV.2 shows cases HMG1, HMG2 and HMG6 which differ 
emily in cross range X, — Xx) and line of sight direction cosine 
value UL but have the same total range. In these cases only 
the far-field beam pattern can cause differences in the magni- 
tude of the frequency response. The table shows that as the 
cross range increases the magnitude falls off, which is con- 
Sistent with the fact that when a beam pattern is steered 


towards end-fire (to larger values of UL, and V,,) the 


Oo 


directivity index decreases due to the increasing beamwidth 
of the main lobe. 

Table IV.3 shows the difference in magnitude of the fre- 
quency response at the center of the receive array between 
HMG1, HMG5, and HMG7. Cases with the same geometry but a 
different value for the single frequency component. This 
difference in magnitude is due to the frequency dependent 
amplitude factor inherent in the WKB approximation, which ds 
not exact for a homogeneous medium. 

Table IV.4 compares the magnitudes of the frequency 
response between the homogeneous cases HMG1, HMG2, and HMG5 
and their inhomogeneous counterparts INHMG1, INHMG2 and INHMG5. 
We see that the magnitudes for the inhomogeneous cases are 
slightly lower, an expected phenomena because of the channel- 
ing effect in the inhomogeneous medium. In a medium with a 
constant positive gradient, the acoustic rays are bent upwards 
and the sound intensity decreases with increasing depth. 

This causes a lower magnitude of the acoustic signal then 
expected by spherical spreading only. 

We will now present numerical and graphical output from 
the three dimensional DFT beamforming program. The frequency 
transform for the single frequency component cases simply 
reveals the magnitude of the transfer function multiplied by 
the amplitude of the input electrical signal and is not shown 
here. Figures 7 through 10 show the magnitude of the spatial 


DET W.2 JG, oer You (G05 ,n) versus direction cosine u, and 
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the magnitude of the spatial DFT w.r.t. y, li.e., icin, So 


“SN 
versus direction cosine u, for the cases HMG1, HMG2, INHMG1L 
and INHMG2. These graphs allow us to estimate the direction 
Of propagation of the incoming plane wave by finding the 
coordinates of the maxima on the graphs. Tables IV.5 to 

IV.8 provide numerical values for the above mentioned DFT's 
w.r.t. x and y. They show values of You (4,0 ,-n) and You (G-m,s) 
around the maxima and make it possible to form more accurate 
estimates. 

From the estimated values uw. and Vo, we can calculate 
estimates of target elevation angle 8 and target bearing 
angle Vo: Table IV.9 gives the results for the different 
cases. The values of the estimators a and v. were determined 
by performing a second order interpolation on the data given 
in tables IV.5 through IV.8. It is seen that the estimates 
@@ectarget angles for the homogeneous cases are correct and 
Saoletne line of Sight angles from the receiver to the trans- 
mitter. The estimated target angles for the deterministic 
inhomogeneous cases deviate slightly from the line of sight 
angles. This is due to the positive sound speed gradient 
which causes upward bending of the acoustic rays so the direc- 
tion cosine Vis directly related to the vertical angle of 
propagation of the plane wave, decreases as the wave travels 
along. The net effect for these cases is relatively small 
Due noticeable. 


The effect of randomness, characterized by a random component 


Mune neomrnGgex GL retraction with a standard deviation of 
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go = 1.0 x1077 is examined next. Figure ll shows estimates of 
ue and Ne for the case RINHMG1 and Fig. 12 shows estimates of 
ve for cases RINHMG1 and RINHMG2. As can be seen, the esti- 
mate of Us is not influenced, which follows readily from the 
fact that the randomness is only present in the y direction. 
The estimate of Vv, as shown in Fig. 12 for the case RINHMG2 
shows the same random fluctuations as the one shown for case 
RINHMG1L in Fig. ll, since both cases used the same seed for 
the random number generator, resulting in the same random 
number sequence. Figure 12 also shows RINHMG1 with a different 
seed for the random number sequence. Although both are only 
two possible samples they clearly show the strong effect of 
the randomness on the beamformer output and the deterioration 
of the guality of the estimates (see Table IV.9). 

Two three-dimensional plots are shown in Fig. 13, one for 
the deterministic inhomogeneous case INHMG2 and one for the 
random inhomogeneous ease RINHMG2. These graphs illustrate 
the kind of output that can be obtained from the DFT beam- 
former, although 1t 1s difimeult “ee aceeduarel ee acmees 
numerical values. 

Finally, the cases HMG8, INHMG8 and RINGHM8, which involve 
several frequency components, are considered. The output 
for these cases is given in numerical and graphical form. 
Table IV.10 compares the Fourier series coefficients of the 
output electrical signal with the Fourier series coefficients 
of the transmitted electrical signal» Wihe widelywevarying Value. 


are caused by: 
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mee ine wbeam pattern Of "tme Eransmit array is frequency 
dependent. The differences between the transmitted frequency 
components are relatively large and cause large differences in 
amplitude of the Fourier series coefficients. The maximum 
values are observed at f = 3000 Hz, 1.e., the frequency for 
which beam steering was done. 

2) The deterministic inhomogeneous medium: As already 
shown in other test cases the amplitude of the received signal 
is frequency dependent due to the WKB approximation to the 
wave equation in an inhomogeneous medium. Eguation (2.10) 


clearly shows the frequency dependence of the medium transfer 


monet lon. 
3) The random inhomogeneous behavior of the medium: The 
randomness introduces additional deviations. teecan. Oo. seen 


from table IV.10 that the effect of the randomness increases 
with increasing frequency. Also, in the case of a random 
medium, the estimation of the Fourier series coefficients 
becomes strongly dependent on the location of the element in 
the receive array. Note that table IV.10 only shows values 
for the receive element (m = 0, n = QO). 

The received angular spectrum for the different frequency 
components of the cases INHMG8 and RINHMG8 is shown in Figs. 
mieehrough 19. For the deterministic inhomogeneous case, 
Bags. 14 through 16 show a regular behavior of the angular 
spectrum according to the rectangular amplitude weighting 


used for the receive array. Estimation of Un 1s exact. 
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Estimate leume: Vv, returns a Slightly different value from the 
line of sight value for direction cosine i due to the 
inhomogeneous behavior of the medium. Figures 14, 15 and 
16 demonstrate the increasing directivity of the receive 
array with increasing frequency, thus leading to more accurate 
estimates for higher frequencies. 

Figures 17, 18 and 19 clearly show the strong influence 
of the randomness on the resulting angular spectrum. The 
randomness creates strong "false" sidelobes and causes an 
offset of the estimated value of direction cosine Ve The 
effects of the randomness increase with increasing frequency. 
At f = 5000 Hz (Fig. 19) the estimate of ve is determined by 
a spurious side lobe and returns a totally wrong value. 
Table IV.11 shows estimated values of Uy and Vyas at the dittiew. 
ent frequencies and their associated estimates of target 


nN 


elevation and bearing angles, OG and Vor respectively. 
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PAbibs fy. 1. 


Magnitude of Frequency Response Comparing Cases 
Which Are Different Only in Total Range 


Case Total Range [shoe Ol (0) 

_ ee (in m) 5 

HMG1 2000 0.02434 

HMG3 1000 0.04868 

HMG4 SUOo OROne: 2:2 
AB hi eed ee cz 


Magnitude of Frequency Response Comparing Cases 
Which Are Different Only in Cross Range 


Case Total Range Cross Range pe (2 0) 0) | 
= (Cabot any), (in m) 

HMG1 ZVI O 0.0 0.02434 
HMG2 Z0010 Z00%0 0202 41:7 
HMG6 2000 eZ 24 7 Oe Oo s4 
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TABLES Ly .3 


Magnitude of Frequency Response Comparing Cases 
Which Are Different Only in Frequency 





Case Frequency i, Once 
(Can SHz) 
HMGL 1000 0.02434 
HMG5 2000 0.03441 
HMG7 4000 0.04866 
TABLE IV e4 


Magnitude of Frequency Response Comparing 
Homogeneous and Inhomogeneous Cases 


Case Hee oO) 

HMG1L - INHMG1 0.02434 - 0.02341 
HMG2 -— INHMG2 0.02417 - 0.02326 
HMG5 - INHMG5 0. O34418=020353 1c 
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ees ey ales S 


MaqyreGiicloe "@mmmenc ya) VS. Direction 
Cosine u for HMGlL, INHMGL 


Direction HMG1 INHMG1 
Cos tae uy os ae 
or 15 2OP2 L416 -0116804 
[0 enks, Os Z eS aro OG o> 
=0'.005 POR Gao Oe -0 3.0 
ie) nO 2 EG yO ILE Ste 
v-6:05 Ow 2 Go. OMe 38 
oe O10 CEZAR SY Sis: ~OubiG 9 45 
O015 One eis .0116804 


ee oe ee 


Magmaeuaenwor YoulG,r,n) VS. Direction 
Cosine u for HMG2, INHMG2 


Direction HMG2 INHMG2 
Cosine 7 <4 

Wer uisS POZO S73 POS OC 
02-090 .0120744 sO ol4 7 
0, 0.95 sOUAO iG SO eG 2S 
0-106 BOZO eige nO G27 4 
oe os SO2 G3 42 .0116254 
Oe e1O OZ 07 516 20 LLG 
Cee 5 20. 20 Gale -0116046 


oo 


TAB IGE sky 5 7 


Magnitude of Ysniq,m,s) Vs. Direction 
Cosine v for HMG1l, HMG2 





WF 


Direction HMGL HMG2 
Cosliew. oe 
0.484 Oe Ge -0120534 
0.489 ~O PZ isan 2OL2Z069'6 
0.494 2002636 .0120800 
0.499 sO LZ esr .0120844 
0.504 O12 leoe POU Oe 
02509 Olga So .0120754 
0.514 .0121457 SO Z0GZ8 
TAB ii ly ae 
Magnitude of Yon (q,m,s) Vs. Di Gee eon 
Cosine v for INHMG1L, INHMG2 

Direction INHMG1 INHMG2 
Cosine v 
0.475 OLb6 jo8 .0115947 
0.480 sOLL6E893 a0) sy aa 
Oxv4iss -O11L7001 Od Zies 
0.490 O17 0S .0116264 
0.495 -OLL7045 Ona Zao 
Oca 0 0 0 Lies s3  OTEGEGS 
0.505 .0116864 -O1Z6073 








eee ies al .<9 


Estimated Values Uo? Vor oe and es 
_ “0 nO Gan fas na ee aaues 
HMG1 Oe 000 Os 500 +2070 ooo 
HMG2 Oo 0500 ono S +10. 02 
INHMG1L 0.000 0.492 +29 .47 Poo eoy 
INHMG2 Ore Oe: e497 te rele +78.40 
RINHMG2 Ore On: 02525 oe ee +79.07 (seed #1) 
RINHMG1L Ors000 Omi S +31.53 +90.00 (seed #1) 
RINHMG1 0.000 Or 0's oor +90.00 (seed #2) 
1) The line of sight angles from receive to transmit 


array (Or 


HMGL, INHMGL, RINHMG1: 9 


HMG2, INHMG2, RINHMG2: 0 


2) Note that 
estimates 
U6 and Vo 
eled ola ty 


target angles) are: 


2 0nO0 sane) +201-00° 


ag 


= ooo 


+0105, ve 


equation (4.14) which is used to calculate 
of target bearing angle from estimates of 
returns two results. This is due to the 
of a planar array. For the above table the 


appropriate values were selected. 


re 


TABLE iV .4¢0 


Estimation of Fourier Series Coefficients 
INHMG8 and RINHMG8 








Yon (G,0,0) for HMG8, 
Pregw@e me, HMG8 INHMG8 
(in HZ) 
1000 O20CTO a7 0.0009705 
2000 O-009s81zZ6 0.0098901 
3000 0.0209 3i66 0.02 CiS 06 
4000 0 01S 8947 0 Ol 220 52 
5000 020024525 QV: 0083476 
Note: The Fourier series coefficients 
mitted electrical signal are Cy 
ke Mh 2s ae 
TABLE i.e 
Estimated Values of u,v_, 8 
Oo O O 
Case Prequeney, Uy ee 
Cinch 2) (in 
HMG8 1000 0.10 Oe 
: 2000 OG 0 250 
se 3000 Oke O50 
" 4000 O21 0.50 
. 5000 OIL: 0.50 
INHMG8 TOC0 O40 0.49 
2000 07.6 0.49 
e 3000 Oe O29 
sf 4000 O21 0.49 
e 5000 0.10 0.49 
RINHMG8 1000 Oke O65 
i 2000 Ore O29 
3000 Oke 053 
4000 Oe Oro 2 
: DOUG O20 7-00-43 
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V. CONCLUSIONS AND RECOMMENDATIONS 


Computer simulation of a mathematical model of wave propa- 
gation through a random, space-variant, time-invariant, ocean 
medium between two planar arrays has been accomplished. The 
computer simulation incorporates an index of refraction which 
is only a function of depth. The deterministic component of 
the index of refraction represents a sound speed profile with 
constant gradient. The random component is assumed to have 
statistical properties independent of depth. The effects of 
the randomness on the electrical output signals at the differ- 
ent receiver elements are assumed to be uncorrelated. 

The computer program was written in FORTRAN. It implements 
both analytical and numerical integration techniques. The 
results from both techniques are in close agreement. 

A number of test cases have been run. Interpretation of 
the data for simple cases show results which are consistent 
with expectations from theory. The model correctly demonstrated 
the influences of the far field beam pattern of the transmit 
array, the source-receiver geometry, and the deterministic 
inhomogeneous behavior of the medium. Also, a random inhomo- 
geneous medium was shown to affect the wave propagation and 
distort the frequency and angular spectra, leading to errors 
in the estimation of the Fourier series coefficients, target 
bearing angle, and target elevation angle. A final test case 


RINHMG8, with a transmit signal containing several frequency 


8 8 


components demonstrated a complex dependence on the above 


mentioned factors, making interpretation of the results diffi- 
cult. The effect of the randomness was very pronounced in 
this case. 

The computer program could be further improved by imple- 
menting a more sophisticated method to simulate the randomness 
of the medium. The assumption that the random phase terms at 
the different receiver elements are uncorrelated breaks down 
for small interelement spacings. The computer program can 
also be easily modified to accommodate a more complex sound 
speed profile in which the gradient is a function of depth, 
fOr example. 

Future application of the model could be devoted to 
investigating the propagation of more complex transmit signals. 
One could also investigate the effects of the randomness of 
the ocean medium on signal processing using different sizes 


of receive arrays, 1.e., increasing the number of elements. 
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